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ABSTRACT 

We present the first three-dimensional calculations of Radiation Driven Implosion of 
Molecular Clouds including self-gravity and ionization. We discuss the effects of initial 
density perturbations on the dynamics of ionizing globules and show that the onset 
of gravitational collapse can be significantly delayed for a multiple of the implosion 
timescale. We demonstrate that Radiation Driven Implosion could be an efficient pro- 
cess for injecting disordered kinetic energy into molecular clouds in the vicinity of 
massive stars. 
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1 INTRODUCTION 

Massive stars have a strong impact on the evolution of their 
parental molecular clouds (MCs). Their fast winds and their 
ionizing radiation cause an important feedback of energy and 
momentum into the circumstellar environment. 

One of the most interesting effects massi ve stars could 
have on M Cs is "sequential star formation" (Elmegreen & 
Lada 1977). When massive stars heat their environment to 
equilibrium temperatures ~ 10 000 K via their ionizing UV 
radiation, the ionized region is divided from the surrounding 
molecular gas by a very thin transition layer called 'ioniza- 
tion front'. Due to the high pressure in the hot, ionized gas, 
a shock front is driven into the cooler molecular gas which 
sweeps up a thin, dense layer of shocked material between 
it and the ionization front. This layer is likely to be gravi- 
tationally unstable and thus may form a new generation of 
stars on timescales of several million years. 

Molecular clouds are highly irregular, with observed 
structures from large scales down to even the smallest ob- 
servable scales. There is ongoing debate whether the mor- 
phology of molecular clouds can be described by a frac- 
tal. Broad linewidths in molecular clouds also suggest ir- 
regular supersonic motion, which indicates the presence of 
turbulence. In the presence of ionizing radiation, theoret- 
ical models predict that density en hancements will likely 
be excavat ed by the ionization front (Elmogrcon, Kimura & 
Tosa 199^ ). They will subsequently protrude into the ionized 
gas, and are directly exposed to the ionizing radiation. 

The dynamical behaviour (Radiation Driven Implosion, 
RDI) in such exposed globule s was studi ed numerically by 
Sandford, Whitaker & Klein (^983 |l984|), who pointed out 



the time evolution can be subdivided into two stages: the 
implosion phase, lasting for about 100 kyrs, and, if gravi- 
tationally stable, a subsequent quiet, quasi-static cometary 
phase, in which the globule is slowly eroded by ionization on 
timescale s of se veral Myrs. Through semi- analy tical models, 
Bertoldi ( |198E| ) and Bertoldi & McKee ( |l990| ) constrained 
the ranges for mass and ionizing flux in which gravitational 
instability in the compressed globules should be expected. 
More recently, HoUenbach & Bally (1998) and Lopez-Martin 
et al. (2001) investigated the photoevaporation of dense 
clumps of gas by an external source. The hydrodynamics 
of photoionized columns in the Eagle Nebula has been mod- 
eled by Williams et al. (2001). 



Sugitani, Tamura & Ogura ( 1999 2000) observed a sam- 
ple of 89 bright-rimmed clouds (BRCs), apparently objects 
with properties as mentioned above, and found an excess in 
the luminosities and luminosity-to-cloud ratios of embedded 
IRAS sources when compared to sources in isolated globules, 
which is an indication for enhanced star formation in their 
BRCs. Near-IR imaging of 44 BRCs revealed that young 
stellar objects seem to be aligned along the axis towards the 
ionizing cluster. There is even evidence for an age gradient 
with older stellar objects closer to the OB cluster and the 
youngest objects well inside the globules, aligned with the 
IRAS sources. This result is consistent with sequential star 
formation while the shock front advances further and further 
into the molecular cloud. 

The interplay between dynamics and self-gravity in ex- 
posed globules has not yet been treated in detail. Using 
the new, fast SPHI method described in Kessel & Burk- 



ert (200C), we present 3D calculations of imploding clumps 



that the mass distribution of stars formed in a perturbed 
molecular cloud affected by RDI may sensitively depend on 
the spatial scales of the perturbatio ns. A similar approach 
was followed by Lefloch & Lazareff (1994). In their models. 



including self-gravity. We find evidence that initial small- 
scale density perturbations can lead to a substantial stabi- 
lization against gravitational collapse during the implosion 
phase and the subsequent cometary phase, even if the pa- 
rameters of the globules suggest them to be highly unstable 
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accor ding to the analytical estimations of Bertoldi & Mc- 
Kee (|199C|). 



2 INITIAL CONDITIONS 
2.1 Density distribution 

The initial conditions are static solutions to the equa- 



tion of Em den (L907), called "Bonnor-Ebert spheres" (Bon- 



Dor 1956). In these spheres, which are supported by an ex- 
ternal pressure at the truncation radius rt , the gravitational 
force is everywhere balanced by the pressure gradient, which 
is vanishing at the center. The scale-free form of the differ- 
ential equation describing this situation is 



: exp(-?/>). 



(1) 



The relation between the parameters ^ and tp on the one 
hand and physical quantities density p and radius r on the 
other hand is given by 

p = Aexp(~^), r = l3'^^\-'^^^, 

where A is the maximum density for tp = and 

kT 



is the temperature of the isothermal gas. Here, k is Boltz- 
mann's Constant, G the Gravitational Constant and fi the 
mean molecular weight. 

Stability analysis shows that there is a critical value 
^pc for the truncation radius of the sphere. Below ipc, the 
sphere is stable against radial density perturbations. Above, 
it commences gravitational collapse on a few gravitational 
timescales if perturbed. At ipc, the sphere contains about 
one Jeans mass of matter, and the density contrast from 
the center to the pressure boundary is ~ 14. Since we don't 
expect such supercritical cores to survive long enough in 
molecular clouds to be influenced by ionizing radiation in 
relevant number, we instead restrict ourselves to the critical 
case with the truncation radius t/jt = tpc- This enables us to 
study whether compression by ionizing radiation can force 
one initial Jeans mass into gravitational collapse. If not, we 
can conclude that induced collapse for cases with masses 
less than a Jeans mass should also be highly improbable. 
Additional simulations of initially Jeans-stable clumps with 
tpt < tpc will be studied in a subsequent paper. 

The maximum particle density of H atoms in the cores 
is assumed to be A = 1000 cm~'^ and the equilibrium tem- 
perature in the neutral cloud material Tc = 10 K to mimic 
the typical conditions observed in molecular clouds. 

The integration volume has the shape of a cube, which 
contains the Bonnor-Ebert sphere, embedded in a surround- 
ing gas of constant density Us- This gas provides the bound- 
ary pressure for the sphere. The diameter of the sphere is 0.8 
times the size of the box. Note that this configuration alone 
would be unstable against gravitational collapse due to the 
additional mass of the surrounding gas, in which pressure 
gradient and gravitational force are not in balance. How- 
ever, the timescale of this collapse is much longer than the 
dynamical timescale which is imposed on the system by the 
large sound velocity of about 13 kms~^ in the hot, ionized 
medium. The setup chosen is not thought to provide a stable 



starting configuration. Our intention is rather to start with 
a reasonable density distribution, in which the initial core 
contains about one Jeans mass. 

We allow material to stream freely away from the initial 
integration volume, i.e. no boundary conditions are imposed. 



2.2 Initial density perturbations 

To investigate the effects of density irregularities, we impose 
in some calculations Gaussian random perturbations on the 
initi al den sity distribution as described by Klessen & Burk- 
ert ( 200C ) . The density distribution is determined by the 
power spect rum P (k) oc k" . This is achieved by applying the 



Zel'Dovich (1970) approximation. Besides the unperturbed 



case without and with self-gravity (in the following cases A 
and B), we ran perturbed cases with a — —3 (more power 
on large scales, case C) and a = — 1 (more power on small 
scales, case D). The amplitude of the perturbations was cho- 
sen to be highly non-linear, with a mean relative variance of 
(Jr = a/p = 0.5. 



2.3 Parameters of the radiation field 

The ionizing flux falls perpendicular onto one of the cube 
surfaces with a rate of photons per time and area, _Fo. We 
specify the "strength" C of the radiation field as the ratio 
between the cube size and the path length I which the radi- 
ation can penetrate into the surrounding gas via an r-type 
ionization front. This length is given by the balance between 
Fo and the recombinations along I: 



n(s) aBds = n^aBl, 



(2) 



where as denotes the recombina tion coefficient adopting the 
"on the spot" approximation (Baker & Menzel 1962). In 



other words, the parameter C determines how much of the 
density enhancement is excavated by the ionization front 
before it turns from weak r-type into weak d-type. From that 
moment on the medium dynamically reacts to the pressure 
difference across the ionization front. For the calculations in 
this paper, we chose C — .97. 



3 RESULTS 

3.1 Time evolution without self-gravity: case A 

The time evolution in case A ( Figs. is analogous to the 
models of Lefioch & Lazareff (1994). One can distinguish 
between four different phases: 



(i) Immediately after the ionizing fiux is switched on, a 
weak R-type ionization front is driven supersonically into 
the medium. After ~ 1.5 kyrs, the denser part of the cloud 
is excavated and surrounded by a region of ionized gas. The 
ionization front stalls and becomes nearly sonic. It converts 
into a weak D-type front. The gas can now react dynamically 
to thepressure imbalance between ionized and unionized gas 
(Fig. |l|, first panels). 

(ii) In the second phase, this imbalance drives a shock 
front into the cloud with velocities ~ 5km/s. The shock 
accumulates a thin layer of dense shocked material (Fig. ^ 
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Figure 1. Model A (see figl.jpg). Radiation impinging from the bottom. Left column: 3D projection. Ionized gas is shown in blue, 
neutral gas in red. Brightness is a measure for the column density. Right column: cuts through the model along the symmetry axis. Color 
coded: log density in g/cm"^. Contour lines: ionization fraction. Arrows: velocity field. h (*'m,i) corresponds to longest black (bright) 
arrow. Evolution time indicated in the upper right corner. 



Figure 2. Model A (continued, see fig2.jpg). Radiation impinging from the bottom. Left column: 3D projection. Ionized gas is shown in 
blue, neutral gas in red. Brightness is a measure for the column density. Right column: cuts through the model along the symmetry axis. 
Color coded: log density in g/cm^. Contour lines: ionization fraction. Arrows: velocity field. h {um l) corresponds to longest black 
(bright) arrow. Evolution time indicated in the upper right corner. 



second panels). It thereby converges towards the axis of the 
finger-like structure ( "globule" ) due to the curvature of the 
ionization front. 

(iii) After 130 Myrs, the globule reaches a stage of maxi- 
mum compression (Fig. |l|, third panels). Caused by the over- 
shoot in pressure, the cloud starts reexpanding and 

(iv) enters a fourth phase, in which gas is being ion- 
ized continuously away from the surface facing the ionizing 
source (Fig. ^). This leads to a mass loss of ~ 10~* Mq/jt. 
Gas is leaving the ionization front with velocities of ~ 
13 km/s, causing a rocket effect which accelerates the gas 
in the tip into the direction of photon flow. During this ac- 
celeration, gas is collected inside the globule in a lid-shaped 
feature facing the ionizing source. It is subject to dynami- 
cal instabilities which lead to the disintegration of the cloud 
after ~ 600 kyrs. The remaining fragments are then ionized 
completely after a few 10 kyrs. 



3.2 Time evolution with self-gravity, no initial 
density perturbations: case B 

Placing our models with R^^/Saq ~ 4/71 and ~ 40 into 
the Figure 14 in Bertoldi & McKee (1990) (reprinted here as 
Fig. ^), they lie far in the regime where gravitational insta- 
bility is expected. Simulations with smooth initial conditions 
including self-gravity verify this conclusion. The overall evo- 
lution is indistinguishable to the one in case A, until the 
point of maximum compression is reached. At that moment 
gravitational collapse sets in, which can be seen in the sud- 
den rise in the maximum density of the system (see Fig. «h. 
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3.3 Time evolution with self-gravity, large-scale 
and small-scale initial density perturbations: 
cases C and D 

In these models, the evolution fundamentally differs from 
the one in model B. During the compression phase, den- 
sity enhancements are amplified in the shocked layer, simi- 
lar to wh at is observed in the simulations by Elmegreen et 
al. (|1995|) (Figs, pi 0) and those of Wilhams et al. (2001) 



(Fig. 6f). They are reminiscent of the observed small-scale 
structures seen in M16 (Hester et al. 1996; White et al. 
1999). Contrary to case B, no collapse occurs at the moment 
of maximum compression. We rather observe an evolution 
analogous to the phases 3 and 4 of case A. After a delay 
of 300 kyrs, the layer starts fragmenting by gravitational 
instability. 



Figure 6. Maximum density (upper panel) and ionization rate 
(lower panel) in the systems. One can dist inguish between the 
four different phases mentioned in Sect. [j.l] . The dense shocked 
layer builds up for t < 130 kyrs. At t ~ 130 kyrs, the systems 
reach maximum compression. Here model B suffers gravitational 
collapse, which reveals as the steep rise in density. For models 
A and C, after a short reexpansion phase there follows a fairly 
quiescent phase, in which material is being constantly ionized 
from the surface. This becomes evident in the lower panel. Note 
the solid line in the upper panel, showing the critical maximum 
density which the numerical method can handle in the cases with 
gravity (Bate & Burkert 1997). 



4 DISCUSSION 
4.1 Jeans Criterion 



As already pointed out in Chapter 2.1, the initial conditions 
were chosen such that the integration volume contains ^ 1 
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Figure 3. Model C (see fig3.jpg). Radiation impinging from the bottom. Left column: 3D projection. Ionized gas is shown in blue, 
neutral gas in red. Brightness is a measure for the column density. Right column: cuts through the model along the symmetry axis. Color 
coded: log density in g/cm"^. Contour lines: ionization fraction. Arrows: velocity field, h (*'m,i) corresponds to longest black (bright) 
arrow. Evolution time indicated in the upper right corner. 



Figure 4. Model C (continued, see fig4.jpg). Radiation impinging from the bottom. Left column: 3D projection. Ionized gas is shown in 
blue, neutral gas in red. Brightness is a measure for the column density. Right column: cuts through the model along the symmetry axis. 
Color coded: log density in g/cm^. Contour lines: ionization fraction. Arrows: velocity field. h {'"ml) corresponds to longest black 
(bright) arrow. Evolution time indicated in the upper right corner. 



Jeans mass. Looking at the Jeans mass Afj and Jeans length 



Energies 



and the freefall timescale, 



2n pGm ' 



Tff = 



37r 
32Gp' 



(3) 



(4) 



it follows that Mj as well as Tg drop with increasing 
density. At the moment of highest compression in model 
A, regions with densities higher than a threshold density 
Pm = 5- 10~^'^gcm~'^ contain 5.3 Mq. For this pm, one finds 
Mj = 2.7 Mq, Lj = 7-10"^ pc and Ta = lOkyrs. Compared 
to these values, the compressed region in the tip of the cloud 
contains 2 Jeans masses and is large enough, and Tg is small 
compared to other timescales in the system. One can thus 
expect that the globule is unstable to gravitational collapse 
at this moment. This explains why we observe collapse when 
including gravity in model B. 

One should however use the Jeans masses criterion with 
care, as the cases C and D illustrate. It only considers the 
balance between internal energy and gravity and is thus a 
static approach, neglecting the effects of internal motion. It 
is as such only applicable in very few cases. 



4.2 Reason for different behaviour of case B with 
respective to cases C and D 

Since the initial conditions are chosen such that they contain 
^ 1 Jeans mass for all cases, the stabilization of cases C 
and D cannot be caused by gas pressure alone. The only 
mechanism which differs in the simulations and which can be 
responsible is undirected gas motion, i.e. disordered kinetic 
energy. 

Fig. shows the total amount of kinetic energy iJkin 
relative to the center of mass compared to the gravitational 
energy -Bgrav in the system for case C. During the compres- 
sion phase, more and more material is collected by the post- 
shock layer and accelerated to velocities 5 km/s. As a 
consequence, i?kin grows steadily until the shock converges 
along the axis of the configuration. Then follows a short 
phase of strong dissipation caused by the colliding shocks. 
During reexpansion and acceleration of the globule, there is 
still a residual amount of -Bkin left in the globule. 

The difference between cases C and D with respect to 
model B is that in the latter case the gas flow is strongly 
converging in the tip of the globule, whereas in the former 
the shock front is refracted by the density fluctuations and 
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Figure 7. Energies vs. time in the neutral gas for model C. See 
text for details. 



thus a substantial amount of i?kin is transformed into undi- 
rected motion on scales small enough to support the glob- 
ule. Additionally, the knots formed in the compression layer 
are relatively stable entities, which behave like a swarm of 
bullets. They tend to pass by each other without merging. 
Both effects inhibit gravitational collapse during maximum 
compression for cases C and D. In 2D calculations using 
cylindrical symmetry, these bullets would be represented by 
concentric rings and their cross section for collision would be 
artificially enhanced. It is thus doubtful whether the delay 
could be observed at all in 2D calculations. 

In the late phases in models C and D, all models show 
gravitational instability in the dense layer which was built 
up in the tip of the cloud facing the source, which leads 
to fragmentation and disruption of the cloud (Fig. ^, low- 
est panel). This happens at different evolutionary times for 
several models we calculated, varying between 420 kyrs and 
600 kyrs. We don't show these late stages in Fig. ^, since 
then over 90 per cent of the gas is already ionized. This 
means that numerical resolution is poor due to low particle 
numbers in the neutral gas. A second reason why we doubt 
the reliability of the time when instability occurs is that 
it will strongly depend on the dissipation rate of the undi- 
rected kinetic energy. Since our numerical scheme uses artifi- 
cial viscosity, we expect numerical dissipation which cannot 
represent the real value. However, the important result we 
trust in is the delay in gravitational collapse after maximum 
compression which we observe in all models with perturbed 
initial conditions. 

The delay in the collapse for models C and D was nei- 
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Figure 5. Figure 14 reprinted from Bertoldi &: McKee (L990), which shows the parameter space of stratified cometary globules. Rp 



is the distance from the ionizing source in parsec, ^49 the number of ionizing photons emitted in units of 10'*^ /s, and m' the mass of 
the globule in units of Mq. The solid lines are thresholds of gravitational stability for different temperatures. Globules positioned to the 
right of these lines should be gravitationally unstable. The models presented here are indicated by an asteric. They lie far to the right 
of the T = 10 K stability line. They should thus collapse after stratification. 



ther ob served nor anticipated in earlier work. Bertoldi & 
McKee (1990) derived the conditions for which gravitational 



virial coefficients 



collapse can be expected by analysing the balance between 
the gas pressure in stratified clouds and the external pressure 
and recoil. In their Figure 14, reprinted here as Fig. ^, all 
the models presented here fall far into the unstable region. 
Note that our globules are far from stratified; the timescale 
for stratification is roughly given by the sound crossing time, 



2r/ac 



(5) 



where r is the radius of the globule and Oc the sound velocity 
in the neutral gas. In our case, 

0.5pc 



0.3pc/Myrs 



IMyr, 



(6) 



which is substantially longer than the timescales considered 
here. 



4.3 Properties of the "bullets" in case C 

In Fig. ^, the distribution of thermal and kinetic virial coeffi- 
cients for the six most massive "bullets" in case D are shown 
after 100 kyrs evolutionary time, shortly before maximum 
compression. The coefficients were estimated by taking the 
ratio of internal or kinetic energy, respectively, compared 
to the gravitational energy at a radius, where the bullet is 
by its higher density still well separated from the embed- 
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Figure 8. Kinetic and thermal virial coefficients for the six most 
massive bullets in model C at t=100 kyrs. 



ding medium. The virial coefficients are all well above 10, 
which indicates that the bullets are no gravitationally bound 
entities, but are pressure-confined. This assumption is even 
affirmed due to the fact that all of the bullets are situated in 
the shocked, dense layer which is formed by the overpressure 
in the neighbouring ionized medium and the recoil caused 
by gas traversing the ionization front with high velocity. 
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The hydrodynamic lifetime of the buUets, estimated 
from their extent and the sound velocity in the neutral 
medium, are found to lie between 100 and 400 kyrs, which is 
comparable to the time the shock front needs to penetrate 
the whole globule. 



4.4 Simplifications and caveats of the models 

4.4-1 Equation of state 

In the calculations presented above, we use a very simplified 
equation of state. We assume that the heating and cool- 
ing timescales are short compared to the other important 
timescales, and that the equilibrium temperature is constant 
at Tc = lOK in the neutral and Ti = 10 OOOK in the ionized 
regions. For material in the ionization front, we couple the 
temperature to the ionization fraction by 



r = Tc ■ (1 



Ti-x. 



(7) 



This ansatz is unsatisfactory in several respects. First, the 
temperature of recombination zones is severely underesti- 
mated, since the temperatures of recently recombined gas 
tend to lie in the regime of several thousand degrees Kelvin. 
We expect that this leads to overestimated densities in the 
shadow region "behind" the main clump. We thus expect 
the mass reservoir to be overestimated as well, leading to 
longer lifetimes and increased probability for gravitational 
collapse in the late stages after rebounce. 

It would be desirable to implement a detailed heating 
and cooling model, which takes into account the tempera- 
tures of dust and gas separately. This implies the application 
of a more general method for determining the radiation field 
generated by the massive stars. One possibility is application 
of flux-limited diffus ion, which was proposed by Levermore 
& Pomraning (1981 ) and already su ccessfully implemented 
into grid codes (Yorke & Kaisig 1995). This issue is currently 
being addressed. 



4-4-2 Diffuse recombination field 

The model used for solving the radiation transfer of ion- 
izing photons uses the simplification called "on the spot 
assumption". Lyman continuum photons produced by elec- 
trons falling into the ground state of protons are themselves 
capable of ionizing H atoms. We assume that this happens 
very close to the site of their production and that following 
the exact equations of radiation transfer for the se ph otons 
can be neglected. It was shown by Yorke et al. (1996) that 



this assumption is not valid in special cases where the diffuse 
recombination radiation leads to the penetration of ioniza- 
tion fronts into regions which are shielded from the direct 
illumination by the ionizing stars. Those regions can thus 
be heated to temperatures around 3000 K. As in the para- 
graph above, this again affects the regions where ionized 
gas streams into the shadows and starts recombining (see 
for example Fig. ^ lower panels). 

We expect that the recombination of this material will 
be delayed, if not prevented, by the Lyman continuum pho- 
tons. This can also be addressed as soon as an implementa- 
tion of the Flux-Limited Diffusion Approximation for SPH 
is at hand. 



5 SUMMARY 

In this paper we demonstrate that the physics of Radiation 
Driven Implosion strongly depends on the properties of the 
initial cloud material. In the presence of initial perturba- 
tions, the evolution is highly nonlinear and the fate of the 
imploding cloud is not easy to foresee. We find that gravita- 
tional collapse need not necessarily occur in those objects, 
depending on the presence, size and strength of initial den- 
sity perturbations. The redistribution of kinetic energy from 
the implosion into undirected motion can stabilize the glob- 
ules for a significant amount of time. Unfortunately, it is 
difficult to estimate this delay from the models presented 
here, since it depends on the decay rate of the kinetic en- 
ergy, and the SPH method used here can reproduce this rate 
only insufficiently. 

The delay could be relevant for star formation if the 
globule is stabilized for such a long time that the ionizing 
radiation can totally erode the globule. If gravitational col- 
lapse occurs before, we expect the collapsing regions to de- 
couple from the gas around them. The remaining gas will 
then be further accelerated by the rocket effect from the 
ionization front, whereas the decoupled cores will lag be- 
hind and eventually be overtaken by the ionization front. 
This process could then repeat for several times, as long as 
there is still mass in the reservoir. We expect the outco me to 
be ve ry similar to the observations of Sugitani et al. ( 199t ; 
200C), where young stellar objects are aligned to the axis of 



the globules with an age gradient from the oldest objects far 
away from the globules and the youngest ones still sitting 
inside. By using the numerical technique called "sink parti- 
cles" (Bate 1995), we are confident to be able to reproduce 



this scenario. One of the next steps should be to take into 
account the measured internal motion in molecular clouds. 
We will present calculations with clouds initially supported 
by turbulent motion in a subsequent paper. 
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